Application of Edwards' statistical mechanics to high dimensional jammed sphere 

packings 
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The isostatic jamming limit of Motionless spherical particles from Edwards' statistical mechanics 
[Song et al, Nature (London) 453, 629 (2008)] is generalized to arbitrary dimension d using a 
liquid-state description. The asymptotic high-dimensional behavior of the self-consistent relation is 
obtained by saddle-point evaluation and checked numerically. The resulting random close packing 
density scaling <f> ~ d2~ d is consistent with that of other approaches, such as replica theory and 
density functional theory. The validity of various structural approximations is assessed by comparing 
with three- to six-dimensional isostatic packings obtained from simulations. These numerical results 
support a growing accuracy of the theoretical approach with dimension. The approach could thus 
serve as a starting point to obtain a geometrical understanding of the higher-order correlations 
present in jammed packings. 
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I. INTRODUCTION 



The sphere packing problem in large spatial dimension 
d is related to several important mathematical problems 
in the context of signal digitalization and of error cor- 
recting codes, in particular. It has been investigated in 
detail by the information theory community [1, 2], but in 
spite of such strong interest, the known rigorous bounds 
on packing fractions ip are not very restrictive. For the 
lower bound, the classical Minkowsky result 4> ~ 2~ d [1, 
Chap. 1, Sec. 1.5] can be improved for lattice packings 
4> > 2d2~ d [3], and Ref. [4] discusses a procedure to ac- 
tually construct packings that achieve this bound. For 
the upper bound, Kabatiansky and Levensthcin have ob- 
tained an asymptotic scaling <f> ~ 2~°- 5990 - d [5]. Though 
the 4> values of laminated lattices up to d — 50 seem to 
suggest that there exist lattices where 2 d <fi grows expo- 
nentially with d [1, Chap. 6], it is quite possible for 
this observation to result from pre-asymptotic effects. 
The gap between the known upper and lower bounds 
thus grow exponentially with d. This broad uncertainty 
leaves open the possibility that the densest packings for 
d — > oo may be amorphous. It has indeed been pro- 
posed in Ref. [6] that it could be possible to construct 
amorphous packings that have a density not only expo- 
nentially higher than the Minkowsky lower bound, but 
actually very close to the Kabatiansky-Levensthein up- 
per bound. A better understanding of high-dimensional 
amorphous and lattice packings would help clarify this 
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intriguing issue. 

Dense amorphous packings of hard spheres are pro- 
duced according to a specific dynamical protocol. Typ- 
ically, one starts from an initial random configuration 
of spheres obtained, e.g., by throwing them into a con- 
tainer, then shaking, tapping, or agitating them until a 
jammed structure is obtained [7-17]. In numerical sim- 
ulations, amorphous packings are produced by inflating 
hard particles while avoiding superposition via molecu- 
lar dynamics [18-20], by compressing deformable parti- 
cles [21, 22], or by minimizing the interaction energy of 
soft particles [23-27]. It is an observational fact that 
these procedures, when crystallization is avoided, lead to 
a final packing fraction close to 0.64 in d = 3. This "ran- 
dom close packing density" , which is approximately 10% 
smaller than the density of the best lattice packing in the 
same dimension, is conjectured to be the densest possible 
packing that does not display local crystalline order. An 
important remark for the following discussion is that the 
random close packings are found to be "isostatic", that 
is, their average coordination number z is at the limit of 
mechanical stability z — 2d [19, 20, 24, 25]. A completely 
satisfactory characterization of the amorphous states of a 
system of identical hard spheres is, however, not yet avail- 
able. The definition of amorphous close packed states is 
still matter of debate [24, 28-31], in part because the 
metastability of the jammed amorphous state with re- 
spect to the crystal order leads to a thermodynamic am- 
biguity [32]. 

Classical statistical mechanics provide useful insights 
into the problem of sphere packing in large d, for both 
lattice (see Ref. [33] for an attempt in this direction) and 
amorphous packings [32] . In the limit of large spatial di- 
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mension mean-field theory becomes exact because each 
degree of freedom interacts with a large number of neigh- 
bors [34]. Additionally, because surfaces and volumes 
scale the same way for large d and geometrical frustration 
between the liquid and the crystal persists [35, 36], nu- 
cleation is strongly suppressed. Mctastability effects be- 
come less important, which reduces the definitional ambi- 
guity. It is thus likely for statistical mechanics to provide 
precise information on the behavior of amorphous pack- 
ings in this limit. 

Edwards proposed a volume ensemble statistical ap- 
proach to study the "out-of-equilibrium" nature of 
jammed states [37-42], and a simple mean-field theory 
based on this approach was recently developed [36, 43- 
45]. At the theory's core is the distribution of Voronoi 
volumes associated with particles. Information about the 
packings can be extracted from the probability distribu- 
tion of this "volume function" . The self-consistent inte- 
gral equation for the free volume that results can then be 
solved analytically or numerically to derive a relation be- 
tween the packing fraction and the average coordination 
number. This approach has been used in low dimensions 
to predict the density and other properties of jammed 
packings. 

In this paper, we provide an alternative derivation 
of the probability distribution of the volume function 
that is based on a large d approach. From the the- 
ory, we derive a general relation between the density of 
jammed packings and their average coordination num- 
ber <j) ~ z/2 d , that holds at exponential order in d and 
is consistent with the one found in Ref. [6] using a dif- 
ferent approach. For isostatic random close packings 
we then obtain <f> ~ (4/3) d2~ d , which is denser than 
the Minkowsky lower bound. Comparing the theoretical 
predictions with computer-generated amorphous pack- 
ings shows that, though generally poor, the agreement 
nonetheless increases with dimension. 

The paper is organized as follows. In Sec. II we 
detail the notation and computer simulation approach. 
In Sec. Ill we introduce theoretical method and pro- 
vide a liquid-state derivation of the self-consistent clo- 
sure. The saddle-point approximate solution for the high- 
dimensional limit is given in Sec. IV. Readers uninter- 
ested by the technical details should skip to Subsec. IV B 
and IV C, where the result and its physical interpreta- 
tion are given. In Sec. V we analyze the low-dimensional 
corrections to the result, and discuss the implications for 
d = 3. We conclude by comparing our results with other 
theoretical approaches (Sec. VI) and discuss possible im- 
provements to the theory (Sec. VII). 



II. NOTATION AND SIMULATION DETAILS 

Following standard mathematical notation, we define 
the (d— l)-dimensional volume, i.e., the surface, of the 
(d— l)-dimensional unit sphere at the boundary of the 



rf-dimensional unit ball 

2ir d / 2 

which is related to the volume of the unit ball Vd — 
Sd-i/d. The discussion that follows will consider hard 
spherical particles of radius a = 1/2 and volume V g = 
Vd/2 d , using the particle diameter as unit of length. We 
denote by p — N/V the number density of particles and 
by 4> = pVg the packing fraction. 

Isostatic random packings in d = 3—6 are generated us- 
ing the simulation code of Skoge et al. [20] for N = 500 
particles. The event-driven molecular dynamics simu- 
lations use a modified Lubachevsky-Stillinger algorithm 
to generate jammed hard-sphere packings. The system 
dynamically evolves according to Newtonian mechanics 
until a diverging pressure is obtained. For sufficiently 
large compression rates 7 the compressed fluid falls out 
of equilibrium, which results in a jammed configuration. 
Structural analysis reveals that these jammed configu- 
rations are isostatic and disordered without any sign of 
crystallization [20] . We find that compressing the system 
with 7 = 10~ 3 in reduced units, until the reduced pres- 
sure Z = j3p/p = 10 12 , with the thermal energy 1//3 set 
to unity, is sufficient to reproduce the results reported 
in the original work [20]. Finite size analysis conducted 
for packings in d — 6 further indicates that no significant 
changes are observed for N up to 3000. 

Simulations are also employed for testing the theoret- 
ically predicted distribution functions of contact spheres 
(Sec. V A). For this task, we randomly generate positions 
for z contact balls on the surface of a rf-dimcnsional cen- 
tral ball, and keep only the configurations that present 
no overlap. This approach guarantees that the resulting 
configurations do not depend on the dynamical sequence 
of sphere addition. 



III. THE METHOD 

Before obtaining a high-dimensional form for the 
mean-field theory, we briefly review the volume function 
approach and generalize it to arbitrary d. Note that a 
more complete discussion of the case d = 3 can be found 
in Ref. [43, 46]. 



A. Calculation of the volume function 

Recall that a Voronoi cell is defined as a convex poly- 
gon whose interior consists of the points that are closer 
to a given particle than to any other. We refer to a 
particle j as the Voronoi particle of particle i if i and j 
share a common Voronoi boundary Bij, defined as the 
(d— l)-dimensional surface that bisects the separation 
between the two particles (Fig. 1). For a given direction 
s, the distance k(s) from particle i to the boundary 
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FIG. 1: (Color online) Voronoi construction. The boundaries 
of the Voronoi cell of particle i are illustrated by the thicker 
lines. The Voronoi boundary Bij between the central par- 
ticle i and a Voronoi particle j is a (d— l)-dimensional face 
bisecting the separation Tij. The Voronoi boundary along the 



direction s is given by h(s) 



j j (2 cos ( 



angle between s and fij. The short-dashed (blue) sphere is 
the region f2(c), where c = 2h(s) = / cos 9ij . 



h(s) = min 



r ij' 



'j 



, 2s • fi 



2s ■ fi 



2 cos( 



(2) 



where c = m /cos fty . Operationally, a Voronoi particle 
j is thus one that minimizes rij'/ (2s ■ fij'). Note that c 
is the diameter of a spherical region 

Cl(c) = {(r,e 1 ,9 2 ,...9 d ^ 1 ) I r/cos0i<c}, (3) 

where (r, 6±, 9%, . . . Od-i) are d-dimensional spherical co- 
ordinates and 6\ is the angle between f and s (see Fig. 1). 
If particle j is truly a Voronoi particle, 0(c) should be 
empty given that the central particle i is at the origin. 
The Voronoi volume of particle i can thus be expressed 
as the angular integral 



Wi= ds 



k{s)° 



It follows that the ensemble average 
volume V of a packing 



(V) = 




(4) 

of the total 



(5) 



= NV d (k(s) d ) =N(Wi), 

where the last line results from assuming that the sys- 
tem is isotropic and homogenous, such that the ensemble 
average of particle i and direction s is equivalent to the 
overall ensemble average of the packing. The reduced 
free volume per particle then simplifies to 



w 



(Wi)-V g V d (li(s) d )-V g 



Va 



V a 



(c d ) 1 



(6) 



The key quantity missing in the analysis of the Voronoi 
volumes is the probability distribution for c. We define 



f(c)dc as the probability that li(s) € [c/2, (c + dc)/2] 
for i at the origin. Note that because for hard spheres 
c G [l,oo), by definition, f(c) can only be non-zero over 
this same interval. We also define the inverse cumulative 
distribution function P>(c) that f2(c) is empty of particle 
centers for i at the origin 



P>(c) = l- 



f(c')dc' 



and thus 



/(c) 



We then obtain that 
w = — — I = 



dP>(c) 
dc 



(c d - 1) f(c)dc 



(7) 



(8) 



1) 



dP>(c) 
dc 



dc 



(9) 



c !i - 1 P>(c)dc, 



where the last line is obtained after integrating by parts 
and noting that the boundary terms vanish. Tests of 
this identity on amorphous hard spheres packings from 
simulations, where P>(c) is obtained directly from each 
packing (see Fig. 7 below), confirm the validity of the 
underlying isotropy assumption. 



B. Liquid state derivation of P>(c) 

From Eq. (9), the problem of identifying the packing 
fraction at jamming is transformed into that of identi- 
fying the form of P>(c) from the structure of jammed 
configurations. The language of liquid state theory is 
particularly well-suited for this task because the jammed 
packings have structural features similar to that of high- 
density liquids. 

Consider the TV-particle probability density Pjy(R ) of 
finding the particles 1,2, ... ,N with configuration R N = 
Ri, R.2, . . . , Rat. Unit normalization is set by integrating 
the particle positions over space 



J p N (R N )dn N = i. 



(10) 



The configurational average (or ensemble average) of a 
many-body observable F(R N ) is then 



<P(R iV )> 



J F(R N ] 



Pv(R iV )dR 



,N 



(11) 



The associated reduced n-particle probability density (or 
n-point correlation function) 

oo 

Pn(R") = (<S(Ri-Ru)<5(R2-R 2 )---<5(Rn-R„)> 
= (iv^W/ PAr(Rn ' RAr ^ )dRAr ^ 

(12) 
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is itself normalized to 

p n (K n )dR n = 



m 



(N-n)\ 



(13) 



For a system with translational invariance we can also 
define the n-particle correlation function 



<?7l(Rl2, R-13 • • • Rin) = Pn(R™)/p r ' 

with normalization 

P U 1 J #n(Rl2, R-13 • • • R-ln)<£R-12 • • • Ri„ = 



(14) 



(N-iy. 



(JV-n)!" 
(15) 

The n-particle correlation function reduces to the pair 
correlation function (or radial distribution function) for 
the case n = 2. 

Following the strategy of Refs. [47] and [48] for ex- 
pressing P> in terms of n-particle correlation functions, 



J 



we define 



m(R;ft) = 



0, otherwise, 



(16) 



as a characteristic function of space point R inside an 
arbitrary region SI. Here, we will specialize to the region 
Sl(c) defined by Eq. (3). We also define the characteristic 
function 



N 



J(Ri; Si) = - m(Ri - Ri; SI)] 



(17) 



of all N— 1 particles outside region SI centered at Ri, with 
Ri the center of particle 1. Considering the cumulative 
probability function P>(Ri;ft) as the probability that 
all N — 1 particles are outside of SI gives 



AT f 

P>(Ri;n) = —f^r / J(Ri;0)P A r(R Ar )dR^ 1 

= 1 4t / ™(R2-Ri;ft)p 2 (Ri,R 2 )dR 2 

Pi(Ri) J 

+ , * / m(R2-Ri;0)m(R 3 -Ri;0)p 3 (Ri,R2,R3)rfR2R3 

-^Pl(Kl) J 

AT-l 

= ^(-l) fc F fe (R i; r!), 



(18) 



where the last expression implicitly defines 

'l , fc = , 



F fc (Ri;n) 



sr ^ r)M / (R^ 1 ) n-i m(Ri <; ft)dR 4 
I 



k > 1 . 



(19) 



If the system has translational invariance, then for k > 1 

^ fe fc+i 
F k (Sl) = — J fffe+i(Ri2, • • . ,Ri(fe+i)) J| m(R H ; Sl)dK u 

(20) 

P 2 (ft), F 3 (Sl), . . ., are respectively the probabilities 
of finding a pair, triplet, etc., within SI. Equa- 
tions (18) and (20) show that the problem can be solved 
exactly only if a complete knowledge of the n-particle cor- 
relation function g n (Ri2, R13 • • • Rin) is available to all 
orders. But an accurate theoretical prediction of g n for 
jammed states is still lacking, even for the lowest order 
pair correlation function g 2 . In the following, we use the 
generalized Kirkwood superposition approximation [49] 



and a theoretically conjectured form for the pair corre- 
lation function in high dimensions to simplify Eq. (18). 
These approximations result in a simple factorized form 
of P> (c) that only depends on the packing density p and 
coordination number z. 

For dimensions greater than one, the generalized Kirk- 
wood superposition approximation offers a way to reex- 
press these higher-order correlations in terms of pair cor- 
relations 



<7n(Rl2, R13, • • • , Rin) 



II 92^). (21) 

l<i<j<n 



In order to proceed any further with this analysis, we 
need to approximate g 2 (r) . Following Torquato and Still- 
inger's suggestion [6] and the results of replica theory [32] , 
we postulate that the pair correlation for a jammed con- 
figuration with z average contacts per particle is angu- 
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larly independent and decomposable into contact and where 0(x) is the Heaviside step function, 
bulk contributions as 

gi{r) ~ — 5(r — 1) + 9(r — 1), (22) The superposition approximation then becomes 

pod-l 



ffn(R-12)Rl3,---,R-ln) — J^52(R-li) |] ^(R-j/c) 

2<j<fc<n 



i=2 



i=2 



nsa(Rii) n t 1 + tc— - 1) - ©a - ^-fc)] 

2<j<k<n 1 d 1 



(23) 



E 



2<j<k<n 



P S d - 



-5{R 3k - 1) - ^ 6(1 - i?, fe ) 



2<j<k<n 



The relative importance of the various terms in the bracket comes from the scaling of their integral over the volume 
V of the spherical region f2 with diameter c and noting that V g /V <~ c~ d 



1+ E ^r"7^(-Rjfe - 1) - E + 



2<j<k<n 



,1-1 



2<j<k<n 



(IRl2 ■ ■ ■ (IRi n 



V n-i + yn-2 (n-l)(n-2) z _ T/n _ 2 (n-l)(n-2) ^ 



V 



n— 1 



(24) 



To first-order we can then make the approximation where 

n 

ff„(Rl2, Rl3, • • • , Rln) - J] ffa( R ii). ( 25 ) 



i=2 

which amounts to saying that spheres 2...n are corre- 
lated with the central sphere 1 but not with each other. 
Though crude, this treatment is reasonable in large d, 
because the sphere surface is very large compared to the 
occupied surface and the packing is increasingly ineffi- 
cient. Plugging this result into Eq. (20) gives 



P_ 
k\ 



,9 2 (Rl2)rfRl2 



(26) 



and in the limit N — > oo 

N-l 



^(«)=g(-D*|(/„^)* 

P / 52(r)dr 

Jo, 



fe=0 

exp 



(27) 



Approximating the pair correlation with Eq. (22) finally 
gives a factorized form whose validity should improve 
with increasing dimension 



P> (c) = exp 



-pV*(c)- 



zS*{c) 
Sd-i 



= P B (c)P c (c), 



(28) 



P B (c)=exp[-pV*(c)} 



and 



Pc(c) = exp 



zS* (c) 



(29) 



(30) 



represent the contributions from bulk and contact parti- 
cles respectively, and 



S*{c)= / 6{r- l)Q{c-r/cos6)dr 



S d . 



s(l/c) 



d8(sm6) 



d-2 



(31) 



and 



V*(c) = / 6(c-r/cos0)dr 



5 d 



j-2 /" 



arccos(l/c) 



^(sin6») d - i [(ccos6») 



(32) 



are empty of particle centers (see Fig. 2). Fig. 2 (b) also 
shows that V*(c) is smaller than the volume of f2(c), 
because of the volume exclusion by the central particle. 
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FIG. 2: (Color online) Schematic illustration of (a) S*(c) and 
(b) V*(c). The white sphere is the excluded zone of the cen- 
tral particle. 



IV. LARGE d ANALYTICAL SOLUTION 

Substituting the factorized form Eq. (28) for P>(c) into 
Eq. (9), one obtains a self-consistency relation for w that 
is amenable to further analytical and numerical treat- 
ment in the high-dimensional limit 



= d 



da 



exp 



2 d V*(c) zS*(c) 



(w + l)V d S d - 



(33) 



and note that s*(— x) = s*(— 1) — s*(x). We reexpress 



Vd Sd-i 
S*(c) _ S d - 2 



8*(1/C) 



S*(-l) 

s*(l/c) 



(36) 



Sd-i S d -i" x ' ' s*(-l) ' 
For large d and i>0we can evaluate s*(x) by saddle- 
point approximation. Developing the exponential around 
£ = x and using the fact that the integrand decays rapidly 
in large d allows to extend the upper boundary of inte- 
gration to £ = oo and leaves corrections of order 0(l/d), 



s *(x)= i dt(i-er- 3)/2 

^ e ^iog(i-« 2 ) 

(d-3)a; 2 



2(d-3)(l-^) 



dx 



1 - Erf 



l(d-3)x 2 
2(1 + a; 2 ) 



(37) 



A. Change of variable 

For notational convenience, we define 

/>arccos(l/c) 

s*(l/c) = / d0(sin0) 



d-2 



= f d£(i-e 2 ) (d - 3)/2 (34) 

Jl/c 

and 

/>arccos(l/c) 

v* (1/c) = / d(9(sin (9)<^ 2 [(c cos (9) d - 1] 

Jo 

=^ - S *(l-2/c 2 )] + (35) 

_1 

-s*(l/c), 



2d -2 



C 'l-i 



where Erf(x) is the error function. Note that the last 
expression is only reasonable for x away from zero. 

By using the relations above and noting that w ~ w + 1 
in the high-dimensional limit, Eq. (33) becomes 



w 



dc c d 1 exp 



+ H d (c;w,z) 



(38) 



with 
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H d (c;w,z) = 



1 \c d , 2Ml-l/c 2 r~^ 2 

s*(-l)\w [ 1 ' w(2d-2) 



ws*(-l) 



1 



(d-l)/2 



2d(l-2/c 2 ) 2d -2 



--z )s*(l/c) 
w 



c / zw\ 
d\ 2*) 



(39) 



1.8 
1.7 
1.6 
1.5 

rH 1.3 
1.2 
1.1 
1.0 



10 1 



10 1 

10 :1 



3 5 7 9 11 13 15 



DOB? 



— o— Numerical solution I 
Numerical solution II 
Asymptotic solution 



10' 



10" 



d 



FIG. 3: (Color Online) The values of l/(df w ) obtained from 
(i) the exact numerical solution of Eq. (33) (numerical solu- 
tion I), (ii) the correction to the asymptotic analysis in large d 
obtained by numerically solving Eq. (48) (numerical solution 
II), and (iii) the leading term Eq. (46) in the asymptotic anal- 
ysis (asymptotic solution). The insert shows the dimensional 
scaling of the corresponding volume fraction <f) from Eq. (47). 
The numerical values for d = 3 — 6 are also listed in Table I. 



B. Asymptotic solution 

For analytical convenience we pose that z and the so- 
lution w(z) have a form 



w 
z = 



= fwWQ 

fz Z I 



where f w and f z are polynomial functions of d. 
tion (42) then becomes 



(43) 



Equa- 



Wq 



2d(l-2/iu§) 2d -2 



Wq 

d 



1 



hJw \ 2 J 



(44) 

The necessity to avoid an exponential divergence of the 
last term imposes 



w = 2/zo, 



(45) 



and allows to express f w as a function of f z . In the iso- 
static case z = 2d, i.e., when z$ = 1 and f z {d) = 2d 
simple algebraic manipulations give / t 
lently, 



^, or equiva- 



w = 73 2 
4d 



(46) 



and 



where the last result is obtained by using the simplified 
expression for s*(x) in Eq. (37). A change of variable 
y = c d /w further reduces the expression to 



1=1 dye 

n/w 



-y+H d [(wy)^ d ;w,z] 



(40) 



Because we expect w to exponentially diverge for large d, 
the lower integration limit should rapidly go to zero. For 
any finite y, y x / d — > 1 in high dimension. The equation 
above then becomes 



dy e -v e H^ 1/d ^*) = e ff*(» 1/i ;».») . (41) 
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(47) 



The scaling form of w provides a way to critically 
examine the simplifications done to s*(x). Using the 
rigorous bounds on packing volume fraction for large 
d, 2- d < cj) < 2-°- 5990d [1, Chap 1. Sec. 1.5], and 
w ~ l/tj), we get 2 > w > 2 5990 w 1.5146. Because 



,i/d 



wo, we use the bound on wo to directly check 
that l/w e [0.5,0.66] and 1 - 2/iog G [0.128,0.5]. We 
conclude that the use of Eq. (37) to approximate s*(x) 
in Hd is justified, because its arguments are always pos- 
itive and bounded away from zero. The validity of the 
simplified expression for Hd in Eq. (41) is also verified. 



Assuming Hd is well behaved at the boundaries, the prob- 
lem reduces to finding a solution for the condition 



H d (w 1/d ;w,z) = Q. 



(42) 



C. Finite d corrections 

Assuming that Eq. (42) has only corrections that are 
subdominant, the above analysis also provides the lead- 
ing finite d corrections to w. These corrections can be 
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obtained by numerically solving the condition Hd(c; w, z) 
with the exact expression Eq. (37) for s* in the case 
z = 2d and w = f w 2 d , then computing f w from the 
condition 



H d (2fi/ d ;f w 2 d ,2d)=0. 



(48) 



The results are compared with the self-consistent solu- 
tion of Eq. (33) obtained by direct numerical evaluation 
of the expressions with the exact form for s* in Fig. 3. 
The quality of the match at high d provides a further 
verification of the simplifying approximations made in 
deriving Eq. (42). 

Both treatments agree rather well down tod« 10. The 
change in behavior of the packing fraction then observed 
corresponds to where the decay of P> (c) becomes slower 
than the growth of c d_1 as the dimension increases. This 
inversion results in c ~ 2 maximizing the integrand in 
high dimensions instead of c ~ 1, as in low dimensions. 
Eq. (42) is derived under the assumption that the in- 
tegral is dominated by c ~ 2, which may explain why 
the two curves deviate from each other in that dimen- 
sional regime. Physically, this change corresponds to 
the jammed packings becoming relatively less efficient 
with dimension. In high dimensions, nearby spheres do 
not provide sufficient cover of the sphere surface, which 
results in the inclusion of farther neighbors within the 
Voronoi cell. In low dimensions, the screening is efficient, 
so the integral is dominated by low c values. Because an 
open structure is less prone to particle-particle correla- 
tions, this interpretation is at least consistent with the 
assumptions made in developing the theory. 



V. LOW d ANALYSIS 

The analysis in the previous section is based on the 
approximate high-dimensional expression for P> (c) from 
Eq. (28) . Errors resulting from the superposition ap- 
proximation and the postulated form for g<i (r) should be 
taken into account, in order to improve 4> results in finite 
dimensions. Though a systematic expression for these 
corrections is not trivial to obtain, it is nonetheless pos- 
sible to guess some of their forms by inspection. And by 
comparing their predictions with the properties of pack- 
ings obtained from simulation, we can assess their suc- 
cess. Our low-dimensional analysis expands and general- 
izes the approach of Ref. [43] for the case d = 3 [43, 46]. 



A. Low d corrections of P>(c) 

We first assume that P>(c) can be factorized in bulk 
Pb(c) and contact Pc(c) contributions as in Eq. (28), 
and examine the correlations between the different types 
of particles separately. The contact term has the general 
form 
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FIG. 4: (Color Online) Comparison of the approximations for 
(S*) from Eq. (52) and Eq. (54) with the simulated results. 
The results show that 1/{S*) is proportional to z when z 
is not too large (z < 2d), (inset) The simulation results in 
d — 20 illustrate the linear behavior of 1/{S*) with z. 




FIG. 5: (Color Online) Comparison of the simulation and 
theoretical forms of Pc(c) from Eq. (49) using different (S*): 
(i) Eq. (52) for a low-density (high-dimensional) limit (theory 
I), (ii) Eq. (54) for a van der Waals-like correction (theory II), 
and (hi) direct surface simulations of (S*) (theory III). 



because 



(S*) = / S*{c)f(c)dc= / S*dP c 

11 , Jo (50) 

P C dS\ 



Pc(c) 



(49) 



where the last step is obtained after integrating by parts 
and assuming that the upper limit of the integration (the 
maximum of S*) goes to infinity. The average of the 
available solid angle (S*) can also be calculated directly 
from simulations of surface sphere configurations. In the 
limit of low occupancy, where the particle volume is neg- 
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FIG. 6: (Color Online) Comparison of the simulation and 
theoretical forms of Ps(c) from Eq. (55) using the packing 
fraction from simulation and different (V*): (i) Eq. (56) for 
the low density limit (theory I) and (ii) Eq. (57) for a van der 
Waals-like correction (theory II). 



(b) 4D 




FIG. 7: (Color Online) Evaluation of the factorization approx- 
imation by comparing simulation results for P > in jammed 
systems with the product of the bulk and contact contribu- 
tions also from simulations. 



ligible and correlations are absent, the probability that z 
contact particles lie outside S*(c) is 



Pc(c) 



1 - 



S*(c) 



Sd- 



(S* 



exp 



zS*(c) 
Sd-i 



s, 



d-1 



(51) 



(52) 



This result is the same as Eq. (30) from the previous 
analysis, because an equivalent approximation to the low- 
density approximation was taken in the high-dimensional 
treatment [6]. Surface correlations are, however, partic- 
ularly significant in low d even for low occupancy, and 



therefore the Sd-i volume (sphere surface) should be 
corrected for the non-negligible volume occupied by the 
particles. We define the surface occupied by a sphere at 
contact 



rrocc Q 



7T/6 



(sin( 



\d-2 



dO, 



(53) 



to include a van der Waals-like correction for the occu- 
pied volume 



(S* 



Sd-i 



cocc 
ZD d-l 



(54) 



This expression reduces to the decorrelated form in high 
dimensions, because S2- C i/ Sd-i vanishes faster than z ~ 
d grows. Figure 4 shows that the latter is slightly better 
at low d and that both approximations rapidly converge 
to the (S*) numerically obtained from random distribu- 
tions of surface sphere configurations. 

A similar treatment can be applied to the bulk term, 
which has the form 



Pb(c) 



-V(c)/(V) 



(55) 



with (V*) the average volume available to particles be- 
yond the contact rim. Without taking into account the 
volume excluded by the central sphere, and assuming 
that there are no angular correlation, the large d limit 
is recovered 



(V*) = V/N = 1/p. 



(56) 



A correction a la van der Waals for the correlations that 
arise from the excluded volume by the central particle 
gives a form similar to that of the surface term 



(O 



V-NV g 
N 



V g w. 



(57) 



The expression reduces to the uncorrelated form in high 
dimensions, because the volume of a unit-diameter ball 
grows much slower than the volume per particle at jam- 
ming. 

The contact and bulk corrections above are equivalent 
to those obtained in the original derivation of the theory 
in d = 3 [43]. Comparing the different expressions with 
the simulation results should inform us on the reason- 
ableness of the approximations made in their derivation. 
As can be seen in Fig. 5, even a crude treatment of vol- 
ume exclusion improves the quality of the scaling form of 
the contact contribution [50]. Using the simulated value 
of (S*) does even better. Moreover, the agreement with 
the scaling form steadily improves from three to six di- 
mensions. Figure 6 shows that the bulk contribution is 
also better captured when correlations are included. We 
note however that though the agreement of the bulk scal- 
ing form improves with dimension, convergence is slower 
than for the contact contribution. Higher-dimensional 
comparisons would be useful to confirm the trend, but 
are unfortunately beyond current computational reach. 
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d 


0sim 


0fact 


01owD 


0high D 


3 


0.64(1) 


0.54 


0.64 


0.39 


4 


0.46(1) 


0.39 


0.38 


0.25 


5 


0.31(1) 


0.27 


0.22 


0.16 


6 


0.20(1) 


0.18 


0.13 


0.10 



TABLE I: Jamming packing fraction of frictionless spheres 
at z = 2d. Compression results from various simulation ap- 
proaches 4>sim (see text) are compared with the integration of 
the surface and bulk decomposition of P>(c) (<^fact), the low- 
dimensional approximation (<^i ow d), and the high-dimensional 
approximation (^highD), as introduced in the text. 



Even assuming that we could obtain separate perfect 
bulk and contact scaling forms, the factorization of the 
two contributions itself remains an approximation. The 
corrections to this approximation are not easily directly 
tractable, but the rapidity at which it vanishes can be in- 
directly evaluated by analyzing the structure of jammed 
configurations. Figure 7 indicates that corrections to the 
factorization vanish fairly rapidly with dimension, which 
supports the validity of the high-dimensional treatment. 



• '/'fact: the factorized approximation of -P>(c) = 
Pb(c)Pc{c) is obtained by computing Pb(c) and 
Pc(c) directly from simulations of jammed config- 
urations. 

The first approximation is the crudest, while the last 
one only assumes P> to be factorizable. The latter should 
therefore be the most accurate and become more so with 
increasing dimension. The factorized version </>f ac t indeed 
steadily approaches the simulation results s ; m with in- 
creasing dimension. The agreement of <j) with </> s i m with 
theoretical sophistication from </>highD to 4>i ac t also gen- 
erally improves for a given dimension. Note that the 3D 
result for frictionless spheres of Song et al. [43] (</>i ow d 
in Table I), which agrees very well with the simulation 
value, defies however this last trend. 

Deviations between theory and simulation confirm that 
correlations are not negligible in low dimensions. The 
assumption of the factorizability of P> as well as the 
theoretical forms of Pg and Pc should thus be refined 
by including the contributions from these correlations, in 
order to obtain a systematic quantitative agreement in 
low dimensions. 



B. Low d jammed packing fraction 



VI. DISCUSSION AND RELATION WITH 
OTHER APPROACHES 



Comparing the numerical jamming packing fraction 
obtained under different approximation schemes with the 
direct compression of the system for various dimensions 
provides a final evaluation of the theoretical treatment 
(see Table I). From a simulation perspective, it is reas- 
suring to note that the packing fraction obtained from 
compression simulations agree not only with the previ- 
ously reported 4-6d values [20], but also with the free 
volume estimates from the metastable liquid state [36], 
the zero-shear limit in 4d [51], and the canonical 3d re- 
sults [7, 52, 53]. 

From the theoretical point of view, three different lev- 
els of approximation for P> are selected to solve Eq. (9). 

• ^highD: the high-dimensional approximation is the 
numerical solution to Eq. (33) plotted in Fig. 3 
(numerical solution I). This approximation is only 
valid in the high-dimensional or low-density limit. 
Therefore, <^highD is expected to have significant de- 
viations from the simulation values in low d. 

• 01owd: the low-dimensional corrections 
Eqs. (49), (55) and (57) to P>(c) give 



P> (c) w exp 



2 d V*(c) S*{c) 

wVd ~1¥) 



(58) 



where (S*) is obtained from analyzing simulations 
of surface sphere configurations [54] . This is exactly 
the same approach as what was used in Ref. [43] for 
the 3D case. 



In this section, we discuss the relation of the results 
obtained above with other approaches to the problem 
of sphere packing in large dimension. The present the- 
ory predicts that amorphous isostatic, i.e., random close 
packed, configurations have a unique packing fraction 
4> ~ (4/3) d2~ d . This scaling lies within the known rig- 
orous upper and lower bounds for packings. But because 
these constraints are not so difficult to satisfy, it is more 
instructive to contrast the scaling form with that of ap- 
proaches developed specifically to deal with amorphous 
packings. 

A first set of results for amorphous packings has been 
obtained by analyzing the behavior of a class of simple 
algorithms, such as Ghost Random Sequential Addition 
(GRSA). GRSA is actually able to construct packings 
up to density <j> = 2~ d 7 which provides a nice way of gen- 
erating amorphous configurations up to the Minkowski 
lower density bound [55]. Random Sequential Addition 
(RSA), where one attempts to add a sphere randomly and 
accepts the move only if there are no overlaps, produces 
more efficient packings. But RSA is already too complex 
to be analyzed analytically, so one has to resort to numer- 
ical investigations [56, 57]. The results of RSA are consis- 
tent with <fi = d 2~ d , which is closer to but still lower than 
our scaling form. It is well known that RSA algorithms 
arc not very efficient in low dimensions. A much more 
efficient approach is the Lubachevsky-Stillinger (LS) al- 
gorithm, which is able to construct isostatic packings 
with a density close to random close packing in low di- 
mensions (see Sect. II) [20]. Unfortunately, investigat- 
ing the LS algorithm in the large d limit is not possible. 
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The extrapolation of the scaling form to low dimensions 
4> ~ 2.56 d2~ d [20] is probably not reliable due to the 
change in the nature of packings around d w 10, as dis- 
cussed above. 

In the absence of a direct way to analyze the packing 
algorithms, we resort to an alternate approach that pro- 
vides an estimate of the random close packing density in 
large d. The mean field approach to the glass transition 
known as Random First Order Transition (RFOT) the- 
ory [58] assumes that amorphous packings correspond to 
the infinite pressure limit of long-lived metastable hard 
sphere glasses [32] . The problem of random close pack- 
ing is then reduced to the study of a simpler equilibrium 
problem, that of computing the equation of state of the 
glass. This approach predicts that the system remains 
a liquid [59, 60] up to a certain packing fraction 0dyn, 
where a dynamical transition to a finite pressure glass 
occurs [58]. Inside the glass phase, a huge number of 
glassy states coexist. Applying an infinite pressure on 
these glassy states results in isostatic jammed configura- 
tions with a range of packing fractions, as has been nu- 
merically verified in low dimensions [20, 61]. In general, 
the range extends from a threshold (th) packing fraction 
to the glass close packed (GCP) packing fraction, i.e., 
4> € [0th, 0gcp] [32]. One of the key qualitative results of 
the RFOT treatment therefore differs from the present 
theory's uniqueness prediction. 

Comparing the detailed high-dimensional behavior of 
the various approaches sheds some light on this disagree- 
ment. Several different implementations of the general 
RFOT approach provide scaling forms for its constitu- 
tive densities. 

• Density Functional Theory (DFT) predicts (j>d yn = 
V2^d2- d ~ 4.13d2- d [58]. 

• Mode-Coupling Theory (MCT) predicts (j) dyn ~ 
d 2 2~ d in its full version [62, 63] and (f>d yn = 
2y/2Tied2- d ~ S.2M2- d when using a Gaus- 
sian approximation for the non-ergodic parame- 
ter [58, 62]. It was shown in Ref. [62], however, that 
MCT leads to inconsistencies above <f> <~ d2~ d . The 
scaling dyn ~ d 2 2~ d predicted by the full MCT 
thus seems unreliable. 

• Replica Theory (RT) predicts cj> dyn ~ 4.8d2- d 
for the dynamical transition, and t h ~ 6.26 d 2~ d 
and 0gcp ~ d ln(rf) 2~ d as boundaries for jam- 
ming. 

Interestingly, DFT, RT, the corrected MCT, and 
the present theory all predict <fi ~ d2~ d for the 
glassy /jamming density. It remains a problem of RFOT 
to reconcile the (j) dyn predictions coming from the differ- 
ent approaches (DFT and MCT cannot make any pre- 
diction for the jamming densities). But, for now, all 
RFOT theories give c/) dyn bigger than (4/3)d2~ d , which 
suggest that the system should still be a liquid at the den- 
sity where the current theory predicts it to be jammed. 



If one believes the RFOT results, this discrepancy sug- 
gests that the present theory might still be missing some 
correlations that, although irrelevant to determine the 
overall scaling of the density, are important for the pre- 
cise determination of the prefactor. The replica the- 
ory prediction that jammed packings exist in an inter- 
val (f> g [6.26d2~ d ,d ln(d) 2~ d ], whose width grows with 
dimension might be also encoded in some missing correla- 
tions. Understanding the nature of these correlations, or 
disproving the replica results, could be considerable ad- 
vances in the microscopic comprehension of amorphous 
high- and low-dimensional jammed packings. 

Finally, an interesting byproduct of the current anal- 
ysis is a relation between the volume fraction 4> and the 
average coordination number z in high dimensions from 
Eq. (45) 




This scaling is consistent with the results obtained in 
Refs. [6, 64], which use a rather different approach. 
The relation thus seems well verified (possibly with sub- 
exponential corrections) for amorphous packings and a 
range of known lattice (see [1, 65] for the volume frac- 
tion and the kissing number of a list of all known densest 
packings up to d — 128). 

VII. CONCLUSION 

In this paper, we have obtained an asymptotic high- 
dimensional density scaling of random close packings us- 
ing a statistical theory. The scaling form is consistent 
with that obtained from other theories. Comparisons 
between the numerical simulations and the structural ap- 
proximations support the theory's validity in high dimen- 
sions. 

The theory provides a general method for relating 
the local surface constraint in a jammed packing to its 
global properties. A simple relation between the vol- 
ume fraction and the average coordination number is 
obtained by constraining the contact value of the pair 
distribution function. We note, however, that the cur- 
rent approach is a mean-field theory that neglects un- 
constrained spatial correlations and coordination num- 
ber fluctuations. The latter might allow for better com- 
pactified packing structures and the former may ex- 
ist in an amorphous packing of monodisperse spheres, 
even though the system appears mescoscopically homoge- 
neous. The low-dimensional results suggest indeed that 
the role of spatial correlations can be particularly sig- 
nificant. Higher level coarse-graining, such as explic- 
itly treating the second-layer neighbors, the contributions 
from coordination number fluctuations and spatial cor- 
relations might thus improve the theory's predictions. 

Though the current theory neglects three- and higher- 
body correlations as well as noncontact two-body cor- 
relations, the consistency of our scaling form with that 
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of the other theoretical approaches discussed in Sec. VI 
nonetheless suggests the presence of a certain universality 
in the packings of amorphous spheres in this limit, that 
is, the impact of spatial correlations on jammed pack- 
ings may be greatly suppressed in high dimensions. If 
the high-dimensional asymptotic behavior is indeed re- 
lated to the asymptotic behavior in the low-density limit 
for any finite d, as was proposed in Ref. [6], then di- 
mensional studies could provide further constructive in- 
formation on the nature of the jammed state. To that 
end, the intrinsic inclusion of higher-order correlations in 
MCT, DFT, and RT appears advantageous, but the theo- 
ries' built-in complexity also partly obscures the physical 
origin of these correlations. The dimensional perspec- 
tive provided by the current framework might thus be a 
promising geometrical starting point for identifying the 
nature and quantifying the corrections present in pack- 
ings of any dimension. Our understanding of the experi- 



mentally accessible two- and three-dimensional packings 
would certainly benefit from such an advance. 

Finally, though the scaling form we obtain gives amor- 
phous packings that are denser than the Minkowski lower 
bound, it is still too far from the rigorous crystalline up- 
per bound for bringing any resolution to the lattice vs. 
amorphous packing question. 
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